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I. INTRODUCTION 



In the nonrelativistic QCD (NRQCD) factorization approach, the rate of semi-inclusive 
quarkonium production at large transverse momentum (pr) is given as a sum of products of 
short- distance coefficients and NRQCD matrix elements The short-distance coefficients 
are calculable as perturbative series in the strong-coupling constant a s , while the production 
matrix elements, at least so far, must be determined by comparison with experimental data. 
The field theoretic ingredients that form the basis for the NRQCD factorization approach for 
quarkonium production are the collinear factorization of hard-scattering processes at large pt 
0, IE H| and the decomposition of the factored semi-inclusive production rate into NRQCD 
operator matrix elements and short-distance coefficients The NRQCD factorization 
formalism predicts that the NRQCD matrix elements are universal (process independent), 
and it also leads to a set of rules jjj for the scaling of matrix elements and interactions with 
v, where v is the heavy-quark or heavy-antiquark velocity in the quarkonium rest frame. 
(y 2 ~ 0.3 for the J/ip, and v 2 ~ 0.1 for the T.) The confrontations of these expected 
properties of the NRQCD factorization formulas with experimental data are among the key 
tests of the NRQCD approach. 

The production of quarkonium at large px in pp collisions provides a particularly impor- 
tant test of the NRQCD factorization approach. At a pr of several times the quarkonium 
mass, the quarkonium production cross section is dominated by a process in which a gluon 
fragments into a quarkonium [5|. Other processes are suppressed by at least l/p? in the cross 
section. Furthermore, the process in which the gluon fragments into the quarkonium through 
a color-octet heavy quark-antiquark (QQ) channel dominates because it is enhanced by a 
factor v A /a 2 relative to the color-singlet process 0. It was pointed out by Cho and Wise 0] 
that a quarkonium that is produced through the color-octet process should have a substantial 
transverse polarization. The reason for this transverse polarization is that the fragmenting 
gluon is nearly on its mass shell, and, hence, is nearly completely transversely polarized. 
In the color-octet fragmentation process, that polarization is passed on to a produced QQ 
pair, which then evolves into a quarkonium state. According to the velocity-scaling rules 
of NRQCD pj, the evolution of the QQ pair into a quarkonium state is dominated by non- 
spin-flip interactions, which preserve the transverse polarization. Spin-flip interactions are 
suppressed by at least v 2 . 

The prediction of substantial quarkonium transverse polarization at large pt relies not 
only on the validity of the NRQCD factorization formulas for quarkonium production, but 
also on the universality of the NRQCD matrix elements and the velocity-scaling rules. There- 
fore, experimental measurements of the polarization prediction test many of the essential 
features of the NRQCD factorization formalism. Such measurements are also important 
in that they can discriminate between the NRQCD factorization approach and the color- 
evaporation model, which predicts zero polarization for the produced quarkonium. 

There have been several calculations, based on the NRQCD factorization approach, of 
the polarization of the produced quarkonium as a function of pt- These include calculations 
of J/ip polarization at leading order in a s 0], at next-to-leading order in a s |8], and at 
next-to-leading order in a s , including the effects of feeddown from ijj' and \c states 0; 
calculations of if) 1 polarization at next-to-leading order in a s jl, H, EH ; and calculations 
of T polarization at next-to-leading order in a s [12]. These calculations are generally in 
agreement with each other, although they are subject to large theoretical uncertainties, 
which arise mainly from uncertainties in the color-octet NRQCD matrix elements. The 
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Collider Detector at Fermilab (CDF) data for J /if) and if)' polarization [13] do not support 
either the prediction of substantial transverse polarization at large pt or the prediction of 
increasing transverse polarization with increasing p?. However, it must be said that only the 
points at the highest pt are more than 1.5 standard deviations away from the predictions in 
Ref. 0. In the case of the T polarization, the CDF data 3 agree with the prediction 2 • 
However, the experimental error bars are too large to allow one to draw firm conclusions 
about the presence of substantial transverse polarization. 

In light of the possible discrepancy between the predictions for J /if) and if)' polarization 
and the CDF data, it seems worthwhile to investigate the sizes of uncalculated contributions 
to the theoretical cross sections. Among such contributions are relativistic corrections, 
which begin in relative order v 2 . It is known in the case of quarkonium decay rates that 
such corrections can be of the same size as the leading contribution Investigations 
of relativistic corrections to quarkonium decays of order v 4 fl6] suggest that the order-f 2 
contribution gives the bulk of the correction to the leading contribution. 

In this paper, we calculate the relativistic corrections of relative order v 2 to the short- 
distance coefficients of both the color-octet and color-singlet terms in the fragmentation 
function for a gluon to fragment into a spin-triplet quarkonium state. We carry out the 
calculation at leading order in a s . In computing the short- distance coefficients, we make 
use of the Collins-Soper definition of the fragmentation function in terms of a quantum 



chromodynamic (QCD) operator matrix element [17|. In carrying out this calculation, we 



confirm the results for the short-distance coefficients at leading order in v [Hj]. As a check of 
our methods, we also compute the correction of relative order v 2 to the three-gluon decay rate 
of a spin-triplet quarkonium state and find that it agrees with the results of Refs. (la, [l8j . 
We find that the relativistic correction to the color-octet fragmentation function is large 
and negative, while the correction to the color-singlet fragmentation function is large and 
positive. 

The remainder of this paper is organized as follows. In Sec. |H]we review the Collins-Soper 
definition of the fragmentation function, and in Sec. IIHI we use the NRQCD factorization 
formalism to write the expression for the fragmentation function in terms of NRQCD op- 
erator matrix elements and perturbatively calculable short- distance coefficients. In Sec. I1VI 
we describe the calculation of the short-distance coefficients. In Sees. |V| and EU we com- 
pute, respectively, the short- distance coefficients for the color-octet and color-singlet parts 
of the fragmentation function for a gluon fragmenting into a 3 Si quarkonium state, through 
relative order v 2 . Finally, we discuss the implications of our results in Sec. IVII| 



II. COLLINS-SOPER DEFINITION OF THE FRAGMENTATION FUNCTION 

In this section we make use of the Collins-Soper definition of the fragmentation function 
to write, in terms of NRQCD operator matrix elements and short-distance coefficients, 
the fragmentation function for a gluon to fragment into a quarkonium state. The Collins- 
Soper definition was first used in a calculation of a quarkonium fragmentation function by 
Ma 0. 

Here, and throughout this paper, we use the following light-cone coordinates for a four- 
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vector V: 



V = (V + , V-, Vl) = (V + , V~, V\ V 2 ), (la) 
V + = (V° + V 3 )/V2, (lb) 
V = (V°-V 3 )/V2. (lc) 



The scalar product of two four-vectors V and W is then 

V-W = V + W~ + V~W + - V ± ■ W±. (2) 

The fragmentation function D 9 _+h(z, //) is the probability for a gluon that has been 
produced in a hard-scattering process to decay into a hadron H carrying a fraction z of 
the + component of the gluon's momentum. This function can be defined, in a light- 
cone gaugeMn terms of the matrix element of a bi-local operator involving two gluon field 
strengths [i^ . In Ref . ^1 , Collins and Soper introduced a gauge- invariant definition of the 
gluon fragmentation function: 

d— 3 r+oo 



= 2^(^-1)^-2) / 



x 



(0|G^(0)£t((TU PH (zfc+) o x ) ^(x-) fea G^(0 + ,x-,0 ± )|0) . (3) 



Here, G^ u is the gluon field-strength tensor, k is the momentum of the field-strength tensor, 
and d = 4 — 2e is the number of space-time dimensions. There is an implicit average over 
the color and polarization states of the initial gluon. The parameter fi is the factorization 
scale, which appears implicitly through the dimensional regularization of the operator matrix 
element. The operator Vh(p+,p x ) is a projection onto states that, in the asymptotic future, 
contain a hadron H with momentum P = (P + , P~ = (M 2 + Pj_)/ (2P + ), P±), where M is 
the mass of the hadron: 

V HiP+>0± ) = l#( P+ > °l) + X) (H(P+, ± )+X\. (4) 
x 

The eikonal operator S{x~) is a path-ordered exponential of the gluon field that makes the 
expression Q gauge invariant by connecting the different space-time positions of the gluon 
field strengths: 



£(x ) ba = Pexp 



POO 

+ig / dz~A + (0 + ,z-,0_ 

J x~ 



(5) 



ba 



where g is the QCD coupling constant, and A^i^x) is the gluon field. Both and G^ u are 
SU(3)-matrix valued, with the matrices in the adjoint representation. The manifest gauge 
invariance of expression (jSJ) allows us to make use of the Feynman gauge in order to simplify 
the calculation. 

The definition (J5J) is invariant under boosts along the longitudinal direction ^?|- The 
fragmentation function is defined in a frame in which the transverse momentum of the 
hadron H is vanishes: P± = 0±. However, the fragmentation function is also invariant under 
boosts that change the transverse momentum of the hadron, while leaving the + components 
of all momenta unchanged ^t}- [This can be seen from the fact that the definition © is 
manifestly covariant, except for dependences on the + components of some quantities and 
on the dummy variable x~ .} An explicit construction of such a Lorentz transformation is 
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given in Appendix El However, for the purposes of this calculation, we find it convenient to 
work in the frame in which P± = 0j_. 

In general, the fragmentation function (jSJ) involves the long-distance dynamics of the evo- 
lution of gluon into a quarkonium state and, hence, is a nonperturbative quantity. However, 
we may evaluate the short-distance part of evolution of the gluon into a QQ state H as a 
power series in a s . A convenient set of Feynman rules for the perturbative expansion of 
Eq. is given in Ref. [l2|. For the purposes of this calculation, we need only the standard 
QCD Feynman rules and the special rule for the creation of a gluon by the operator v in 
Eq. The latter rule, in momentum space, is a factor 

+*k + (> - ^) <U (6) 

where k is the momentum of the field-strength tensor, Q is the momentum of the gluon, a 
and a are the vector and color indices, respectively, of the created gluon, and b is the color 
index of the eikonal line. In the absence of interactions with the eikonal line, k = Q. 

In the definition (J3J), the transverse momentum of the hadron H is fixed. Therefore, the 
phase space is the product of the phase-space factors of all of the other final-state particles 
and the light-cone-energy-conserving delta function that is implied by the integration over 
x~ in Eq. (jSJ): 

4ttM / , . ^ ,\ dafd d ~ 2 a ± 

d ^=— s { k+ - p+ -E-tjn^^- m 

where S is the statistical factor for identical particles in the final state, a« is the momentum 
of the zth final-state particle, and the product is over all of the final-state particles except H. 
Since we use nonrelativistic normalization for the hadron H, a factor 2M has been included 
in the phase space in order to cancel the relativistic normalization of H in the definition Q. 
In the remainder of this paper, we use nonrelativistic normalization for the hadron H and 
heavy quark Q and antiquark Q, and we use relativistic normalization for all of the other 
particles. 



III. NRQCD FACTORIZATION 

The fragmentation of a gluon into a heavy-quarkonium state H involves many momentum 
scales, ranging from the factorization scale of the fragmentation function /i, which we assume 
to be of the order of the heavy-quark mass m or greater, to momenta much smaller than 
m, for which nonperturbative effects are large. The NRQCD factorization formalism allows 
one to make a systematic separation of momentum scales of order m and larger from scales 
of order mv or smaller. Following this approach pj, we write the fragmentation function of 
a gluon fragmenting into a heavy quarkonium in the form 

D g ^ H (z,n) = [dn(z,n)(OS) +d' n (z,ii)(V%) }+0(v 3 ), (8) 

n 

where the and are NRQCD operators, and the d n (z,fi) and d' n (z,fj) are short- 
distance coefficients. The first term on the right-hand side of Eq. (JSJ) contains the operator 
matrix elements of leading order in v, and the second term contains the operator matrix 
elements of relative order v 2 . The index n represents the quantum numbers of the operator. 
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The matrix elements (Off) and (Vff) are denned in the rest frame of H and are given, in 
the case of a 3 Si state H, by the vacuum expectation values of the following four-quark 
operators: 



£>f = X VV ^2\H + X)(H + X\^a t 



X- 



x 



o H 

^8 



V? 



x 



VT> ^2\H + X)(H + X\ ^a l T a 



X- 



x 



2m? 
1 

2m? 



xV UdYiP ^\H + X)(H + X\ v»V*x + h. c. 
X VT a (§B) % |# + X) + X| ^ VT°x + H. c. 



(9a) 
(9b) 

(9c) 

(9d) 



where m is the heavy-quark mass, ip is the Pauli field that annihilates a Q, x is a Pauli 
field that creates a Q, and x Dip = x (Dip) — (Dxyip- The sum is over all final states that 
contain the specified quarkonium state H. The NRQCD matrix elements are nonperturba- 
tive in nature, but they are universal, in that the same matrix elements describe inclusive 
production of 3 Si quarkonium states in other high-energy processes. 



IV. SHORT-DISTANCE COEFFICIENTS 

The short- distance coefficients in Eq. (JHJ) are independent of the long-distance dynamics of 
hadronization and, consequently, they are independent of the hadronic state H. Therefore, 
they can be calculated by examining the fragmentation function for the case in which the 
state H is a free QQ state. We take the momenta of the Q and Q to be p = P/2 + q and 
p = P/2 — q, respectively. The heavy quark has three-momentum q in the QQ rest frame, 
and, so, the invariant mass of the QQ state is P 2 = 4E 2 , where E = \J m 2 + q 2 . Since 
the momentum q is fixed, the phase space of the QQ pair is that of a single particle with 
momentum P, and it is omitted in the phase-space integrations in the computation of the 
fragmentation function. The fragmentation function for this free QQ state is 

D^ Q q(z^) = [d n (z^)(0^) + d' n (z,fi)(V^)] + 0(v 3 ) . (10) 

n 

The definitions of the NRQCD matrix elements (0® Q ) and (V® Q ) are same as those in 
Eq. (J2J), except for the replacement H — > QQ. The short-distance coefficients d n (z,fi) and 
d' n (z,fi), which are common to both Eq. (JSJ) and Eq. (jTUjl . can be obtained by comparing 
a perturbative calculation of the fragmentation function on the left-hand side of Eq. (JIUj) 
with a perturbative calculation of the NRQCD matrix elements on the right-hand side of 
Eq. pop. While the fragmentation function and the matrix elements in Eq. (JHJ) may both 
display sensitivity to long-range (infrared) interactions, that sensitivity cancels in the short- 
distance coefficients. _ 

In determining the coefficients of (Off Q ) and (V® Q ) in D g _>Qq(z , /i) , it is convenient to 
make use of projection operators for the spin and color states of the QQ pair. The projection 
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operators for the QQ pair in the color-singlet and color-octet configurations are 

Ai = ^=fy, (Ha) 

As = (Hb) 

where i and j are the quark color indices in the fundamental (triplet) representation, a is 
the octet color index, and N c = 3. Spin-projection operators, accurate to all orders in v, 
have been given in Ref. [16]. 1 For the spin-triplet state, the projection operator is 

P 4- 2E 

A(P, q, e*) = A"(P, q)e* a = N& - m) f^-^^+m), (12) 

where e is the polarization of the heavy-quark pair, and N = [2\/2E(E + m)]" 1 . Note 
that we use nonrelativistic normalization for the heavy-quark spinors. If C is the full QCD 
amplitude for a process, then the spin-triplet part of that amplitude is Ai = Tr(CA). 
The S"-wave part of Ai is 

-Ms™ = -Mo + ^M 2 + 0(q 4 ), (13) 

where the first two terms on the right-hand side of Eq. (fT3*j) are the leading and first sub- 
leading terms in the v expansion. Here, 



2(d - 1) dq a dqP 
where 



Mo = Ai\ q ^ , (14a) 

m 2jaf3 Q2 M 



Ai 



2 



(14b) 



jaP = _ g *f3 + p^pP/^ E ^ . ( 15 ) 

The matrix elements of the QQ NRQCD operators are normalized as 

(O? 3 ) = 2(d-l)N c , (16a) 

(Of®) = (rf-l)(iV c 2 -l), (16b) 

(V^) = 4<^>- (16c) 

The quantities that we calculate in this paper are all finite. Thus, we may put d — 4 in 
Eqs. ©, 0, and (OH). 



V. COLOR-OCTET CONTRIBUTION 

Let us calculate the color-octet contribution to the fragmentation of a gluon into a QQ pair 
in full QCD. The Feynman diagram for the color-octet part of the fragmentation function in 



Projection operators that are accurate to lowest order in the nonrelativistic expansion have been given 
previously |21| . 
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FIG. 1: Feynman diagram for the color-octet contribution at leading order in a s to the fragmen- 
tation function for a gluon fragmenting into a color-octet spin-triplet QQ pair. 

leading order in a s is shown in Fig. ^ The circles represent the the gluon field strengths and 
the double lines represent the eikonal operator. The momentum k = (k + , k~,k±_) flows into 
the circle on the left and out the circle on the right. The vertical line represents the final- 
state cut. In leading order in a s , the final state consists of a QQ pair with total momentum 
P = (zk + , P 2 /(2zk + ), 0_|_). The Feynman diagram in Fig. ^includes relativistic corrections 
to all orders in v. The phase-space element is obtained by setting n = and S = 1 in 
Eq. ©: 

= (17) 

The color-octet contribution can be extracted by tracing over the projection operator Ag 
[Eq. (jllbjl ] on both the left-hand and right-hand sides of the cut. Summing over the color 
index a, we obtain the color factor (AT? — 1)/2. The spin-triplet contribution can be extracted 
by tracing over the projection operator A a [Eq. ([12)1] on the left-hand side of the cut and 
the projection operator on the right-hand side of the cut. Multiplying by the prefactor 
in Eq. Q and the phase-space element (fTTj) . we obtain 



D^{z) 



Ka s MM^ c M P sZ^ 
8E 4 



9va 



9nP 



5(1 -z), (18) 



where J crr is given in Eq. ([T5]). and 



M 



aa 

5-wave 



Tr ( 7 a A c 



' S- wave 



■V2g c 



1 



6m 2 



0(v 4 ). 



(19) 



In the last equality in Eq. (fT9"j) . we have made use of Eq. (JHJ). Contracting vector indices, 
factoring out the normalization of the octet matrix elements given in Eqs. (j!6b|) and ()16c|) . 
and setting M = 2E, we obtain 



7ra s <5(l — z) 
3 (AT? - l)m 3 

d 8 (z) (0^) + d' 8 (z) (V^) + 0(v 



0(v 4 



(20) 



where d$(z) and d' 8 (z) are the short-distance coefficients. We obtain the contribution to the 
fragmentation function for a gluon fragmenting into a quarkonium H by replacing the QQ 
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A 




FIG. 2: One of the Feynman diagrams for the color-singlet contribution at leading order in a s to 
the fragmentation function for a gluon fragmenting into a color-singlet spin-triplet QQ pair. The 
other diagrams are obtained by permuting the connections of the gluons to the heavy-quark lines. 



operator matrix elements in Eq. f!2()l) by H operator matrix elements: 

11 



H na s 5(l-z) H 

^ W " 3(jV 2 _ 1)m3 ^8 ) 



1 - ~-vl + 0(v*) 



(21) 



where 



vl = (K)/(0£). (22) 

Our result for the term of leading order in v in Eq. ()21|) is in agreement with the results of 
Refs. HHJ. 



VI. COLOR-SINGLET CONTRIBUTION 

Now let us calculate the color-singlet part of the fragmentation function for a gluon 
fragmenting into a QQ pair in full QCD. The Feynman diagram for this process at leading 
order in a s is shown in Fig. El We assign the momenta and polarization indices for the 
particles in this process as follows: 

g*(k, a) -> QQ(P, a) + g(a, fi a ) + g(b, fi b ), (23) 

where k, P, a, and b are momenta and a, a, /i a , and /x& are polarization indices. Owing 
to the color-singlet and spin-triplet quantum number of the QQ pair, the final state must 
contain at least two gluons that couple to the heavy-quark line. Hence, in leading order in 
a s , there are no diagrams for this process, in any gauge, in which a gluon couples to the 
eikonal line. 

In the frame in which Eq. (J3J) is defined, we can choose the transverse direction so that 

2 



b 

wk + , - =7 i , b± cos d>, b± sin <b ) , (24c) 
2wk + 1 
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where <p is the azimuthal angle of b relative to a. The variables y and w are the light-cone 
fractions of final-state gluons: y = a + /k + , and w = b + /k + , with z + y + w — 1. The phase- 
space element of the process can be obtained by making the substitutions n = 2, S = 2, 
ai = a, and a 2 = b in Eq. (J7j). It is useful in carrying out the phase-space integration to 
express the variables a±, b±, and <fi in terms of Lorentz-invariant dimensionless variables e a , 
e&, and x: 



P ■ a 1 / za\ y\ 
~P2~ = 2 V(2E) 2 y + I J 



eb = p^ = H(#^ + "J' (25b) 

- = ^ = ^f^ 1 + ^ 1 - 2«^x cos . (25c) 



P 2 2(2E) 2 \ y w 

Note that only the invariant variable x depends on the angle <fi. In terms of the variables e a 
and eb, the phase space is 

M(2E) 4 f , . 

Az 2 k+ 2tt Y V ; 



where 



and we have integrated over the azimuthal angle of a and the variable w. The 9 functions 



in Eq. (|26bJ) impose the requirement of the positivity of the variables a\ and b\. We defer 
the discussion of the change of variables from to x to Appendix IBl 

The color-singlet contribution can be extracted by tracing over the projection operator 
A x [Eq. (jlla|) ] on both sides of the final state cut. Using charge-conjugation symmetry to 
relate diagrams involving permutations of the gluon connections to the heavy-quark lines 
(Fig. EI), we find that the color factor is 

d abc V _ (iV 2 -4)(iV 2 -l^ 



A^W C ) 16iV 2 



The spin-triplet contribution can be extracted by tracing over the projection operator A a 
[Eq. (|12j) ] on the left-hand side of the cut and the projection operator A' 3 on the right-hand 
side of the cut. Multiplying by the prefactor in Eq. (jHJ) and the phase-space element (I26aj) . 
we obtain 

D QQ (z) = W - 4)qgM f Z-A 
1 U 32ttz 5 -<W- 7 $ (l + 2e a + 2e fc + 2x) 2 ' 1 1 

where the scalar product in S ■ A denotes the sum over all the repeated vector indices 
appearing in the product of the tensors E and A. These tensors are defined as 

t ( Kn a \ ( k^n p \ 

- 9 /IV g^LaVa 9^ b U b (TT I 9va 7^7" I I 9^ I , 

A = M®M* = MtZTe MfZ a Z , (29b) 



l a -r^ ^ a ^4 7^A CT (P, q) 

p— (c — m p— p — m 



+ 5 perm. , (29c) 
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where I aT is given in Eq. (JT5j) . and the terms labeled "perm." in Eq. ()29c|) are generated by 
permuting the gluon momenta and polarization indices (—k,a), (a,/i a ), (&,//&). In Eq. (|2*§|). 
we have made use of the gauge invariance of the Collins-Soper form of the fragmentation 
function (jHJ) to carry out the calculation in the Feynman gauge. We obtain the S"-wave 
contributions of leading order in v and of relative order v 2 in Eq. (j29b|) by applying Eqs. (|13j) 
and (fH|) to the spin-triplet amputated amplitude (|29c|) . 

Substituting iV c = 3, using Eq. (JTHJ) to factor out the NRQCD matrix element (Of Q ) = 
6N C , and retaining terms through relative order v 2 in 3 • A and in the phase-space factor 
M = 2E, we obtain 

D QQ (z) = MUL f * + + i c ,„m 

1 1 ; 2592ttz J 9 (1 + 2e a + 2e b + 2xf 1 ; 

= d x (z) {Of 3 ) + d' x {z) (Vf 3 ) + 0(v 4 ), (30) 
where d\(z) and d[(z) are the short-distance coefficients, and 

F = E-(M ®M*)\ q ^ , (31a) 

(31b) 



[Note that, although A4o and M.2 are independent of q, the factors 3 in Eq. (|3*T]) introduce 
a dependence on q 2 .] We obtain the contribution to the fragmentation function for a gluon 
fragmenting into a quarkonium H by replacing the QQ operator matrix elements in Eq. (j3*Uj) 
by H operator matrix elements: 

D „ = 5a>(Of) t F a+ (IF„ + F^l 
lV; 25927T2 (1 + 2e a + 2e b + 2x) 2 y h K J 

where v\ is given by 

(7?>/<0f>, (33) 



rather than by Eq. 

The evaluations of F an d ^2 i n Eqs. ()31|) are straightforward, but quite involved. We 
compute F and F 2 using reduce j24Jj . As a check, we carry out independent calculations 
using the feyncalc package [Hf in mathematica j^. By making the replacements 
3 — > —g^ a v a 9^ b u b 9ai3luT and k 2 — >• and multiplying by the appropriate factor, we also 
obtain the decay rate of the process T — > ggg, including the relativistic correction. Our 
results agree with those in Refs. 0,0. The rather lengthy expressions for Fq and F 2 are 
given in a mathematica notebook file |27| . 

In order to evaluate the integrals in Eq. (J32J), we use Eq. ()25c|) to replace with x, as 
we have already mentioned, and we carry out the y integration analytically. This procedure 
is described in detail in Appendix El The integration methods that we use are similar to 
those in Ref. 0. However, in Ref. j5|], the short-distance coefficient d\(z) is expressed as a 
two-dimensional integral, which is then evaluated numerically. In the present calculation, we 
carry out the integrations over the three variables e a , e&, and x numerically. This procedure 
allows one to extract information about the energy spectrum of the radiated gluons. The 
results of this numerical integration are presented in Fig. El and in Table |1J We note that 
the short- distance coefficient d{(z) becomes negative for some values of z. The complete 
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fragmentation function, to all orders in v, is an integral of the square of a quantity and is, 
therefore, positive. However, the individual contributions in the v expansion need not be 
positive. 

We estimate the fragmentation probability by integrating the fragmentation function 
Di(z) over the longitudinal fraction z: 

/ dz D?{z) = 8.29 x 1(T 4 ■ -§<0f> [1 + 2.45 v 2 + 0{v 4 )} . (34) 
Jo m 

Our result for the term of leading order in v in Eq. (|H4|) is in agreement with the result of 
Ref. @. 



VII. DISCUSSION 

We have computed the contributions of leading order in v and the relativistic corrections 
of relative order v 2 to the fragmentation function for a gluon to fragment into a 3 Si heavy- 
quarkonium state. We have computed both the contribution in which the produced QQ 
pair is in a color-octet state and the contribution in which the QQ pair is in a color-singlet 
state. Our results of leading order in v agree with those of Refs. [22J, |23( for the color-octet 
contribution and with those of Ref. for the color-singlet contribution. Our results for the 
corrections of relative order v 2 are new. 

We estimate the relative sizes of the relativistic corrections for fragmentation into J/if> 
by taking v \ and v\ from the Gremm-Kapustin relation |28( 

v x=^ = — — , (35) 

?TlQCD 



12 



TABLE I: Numerical values of the color-singlet short-distance coefficients d\{z) and d'i(z), which are 
defined in Eq. (|30j) . The scaling factors in this table are c\ = 10~ 4 x a 3 /m 3 and d x = 10 -3 x cv 3 /m 3 . 
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14.7 



where m po i e is the pole mass and wiqcd is the mass that appears in the NRQCD action. The 
Gremm-Kapustin relation follows from the equations of motion of NRQCD and is accurate 
up to corrections of relative order v 2 . In the original work of Gremm and Kapustin j2^|, a 
relation was given only for v\. In Eq. (|35|) . we have included the Gremm-Kapustin relation 
for v$, which can be derived in exactly the same manner as the Gremm-Kapustin relation 
for v\. Dimensional regularization of the matrix elements (Vi) and (Vg) is implicit in the 
Gremm-Kapustin relation [29]. Taking ttiqcd = ^poie = 1.4 GeV and Mj/,^ = 3.097 GeV, 
we obtain vf = v$ = 0.21. One should regard this as only a rough estimate of the sizes of 
v\ and Vs. In fact, for m po i e in the range 1.2 GeV < m po i e < 1.6 GeV, which corresponds to 
the latest Particle Data Group compilation [3(}, the values of v\ and u§ given by Eq. (|3*oT) 
can even become negative. 2 On the other hand, the estimate for v\ that we obtain is in 
accordance with expectations from the NRQCD velocity-scaling rules 0], and it lies in the 
central part of the range 0.03 < v \ < 0.6, which follows from a lattice calculation j3j| under 
the assumption that m c lies in the range 1.2 GeV < m c < 1.6 GeV. Inserting v \ = v f = 0.21 



2 Note that negative values of v\ are allowed since, owing to the subtractions of power divergences that are 
implicit in dimensional regularization, the corresponding matrix element is not positive definite. 
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into Eqs. (j21j) and (|34jh we find that the relativistic corrections change the short- distance 
coefficients for the color-octet contribution and the color-singlet contribution by about —40% 
and 50%, respectively. 

We see that, in the case of the J /if), the estimated relativistic corrections are quite large 
in comparison with the contributions of leading order in v 2 . This is not unexpected, given 
that such large relativistic corrections also appear in charmonium decays. 3 Nevertheless, 
these large relativistic corrections cast some doubt on the validity of the v expansion for 
charmonium. One can hope that, as is the case in charmonium decays 0], the corrections 
of relative order v 4 will turn out to be significantly smaller than the corrections of relative 
order v 2 . 

The value of the color-singlet matrix element {Of) is fixed by quarkonium decay rates. 
Therefore, the large relativistic correction to the short-distance coefficient of the color-singlet 
contribution to the fragmentation function will directly affect the theoretical predictions for 
quarkonium production rates. However, in the case of J/i/j production at the Tevatron, the 
color-singlet fragmentation contribution is less than 5% of the total theoretical prediction 
over a wide range of px 0, Hif . Therefore, the relativistic correction to the color- 
singlet short-distance coefficient will have little effect on the predictions for either the J /if) 
production rate or polarization at the Tevatron. 

The color-octet matrix element (Oq) in Eq. (J2Tj) is, at present, obtained by fitting the 
Tevatron data for J / if> production to the complete theoretical expression, which is dominated 
at the largest values of px by the color-octet fragmentation contribution. Therefore, the 
large decrease in the short- distance coefficient for the color-octet contribution will result in 
a corresponding large increase in the fitted value of the matrix element (0$). The net result 
is that the phenomenology of J / ip production at the Tevatron, for either the pt distribution 
or the polarization, will be largely unaffected by the relativistic correction to the color-octet 
short- distance coefficient. Since, in leading order in a s , the matrix element (0$) typically 
appears only in the fragmentation contribution to a process, we expect the phenomenology of 
other processes to be similarly unaffected. We note that, even with the increase in the fitted 
value of (Og) that would result from the relativistic correction, the value of that matrix 
element would still be somewhat smaller than is expected from the velocity scaling rules. 4 
However, the relativistic correction may be important in comparing phenomenological values 
of (Og ) with future lattice calculations. 
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APPENDIX A: LORENTZ TRANSFORMATION FOR THE FRAGMENTATION 
FUNCTION 



We wish to construct a Lorentz transformation that changes the transverse components 
of the hadron's momentum while leaving the + components of all four-vectors unchanged. 
We will construct a particular transformation from a frame in which the fragmenting gluon 
has vanishing transverse momentum to a frame in which the quarkonium has vanishing 
transverse momentum. Such a transformation was already discussed in Ref. Here, we 
give an explicit construction. 

Given a basis set of four linearly independent four- vectors that span the four- dimensional 
space-time, a transformation is a Lorentz transformation iff it leaves all scalar products of 
the basis four-vectors invariant. Using the notation V = {V + , V~, V 1 , V 2 ), we choose for 
our basis set in the original Lorentz frame 
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p 2 +p 2 ± 



2ztf 
(0,1,0,0), 
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P±,0 



(Ala) 

(Alb) 

(Ale) 
(Aid) 



where we have taken the 1 direction to be along p±. We assume that n is invariant under 
the Lorentz transformation, which implies that the + components of all four-vectors are also 
invariant. We also assume, for simplicity, that e 2 is invariant. Then, the transformed basis 
vectors are 



k' 
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zk~ 



k 2 + p 2 Jz 2 



-P± 
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2k- 
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2zk+ 
(0,1,0,0), 
(0,0,0,1). 



0,0 



(A2a) 

(A2b) 

(A2c) 
(A2d) 



The transformed vectors are completely fixed by the requirement that all of the scalar 
products be invariant under the transformation and the requirements that p' ± = 0±, n' = n, 
and e' 2 = e 2 . 

Using Eqs. (|A1J) and (|A2|) . we can construct an explicit transformation matrix L in the 
light-cone basis n = (1,0,0,0), n = (0,1,0,0), e 1 = (0,0,1,0), and e 2 = (0,0,0,1). L is 
defined by V = LV, where V is any four-vector, expressed as a column vector. It follows 
that L is given in terms of the untransformed and transformed light-cone basis vectors by 
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(A3) 
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2k+ 



From Eq. (jAl|) . we see that the light-cone basis vectors can be written in terms of k, p, n, 
and e 2 as 



(A4a) 
(A4b) 
(A4c) 
(A4d) 
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Then, the transformed light-cone basis vectors are 



-' 1 (v e 



n = n, 



2{p' — zk') + n 



z 2 k 2 — p 2 — p 2 ± 
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(A5a) 
(A5b) 
(A5c) 
(A5d) 



2p 

e' 2 = e 2 . 

Substituting Eqs. (lA"4l) and (1X51) into Eq. (IX3|) and using Eqs. (lA"T1) and (|Al|l . we obtain 

/ 1 0\ 

L = 



2(p+) 2 1 p+ u 



V o 001/ 



(A6) 



The inverse transformation from the primed coordinates to the unprimed coordinates is 
given by L^ 1 = L(pj_ — > — 

APPENDIX B: PHASE SPACE FOR THE COLOR-OCTET CONTRIBUTION 

In this appendix we describe the manipulations that we carry out on the phase space for 
the color-singlet contribution to the fragmentation function [Eq. 1)26)1] in order to put it into 
a form that is suitable for numerical integration. We write the phase-space integration as 



de n 



1-2 



2tt 



(Bl) 



where = (e a - £) (e 6 - £). 

Since the integrand depends on <fr only through the variable x [Eq. ()25cj) ]. it is an even 
function of 0. Therefore, we may restrict the range of to to ir and double the integrand: 
f Q 71 dcj)/(2ir) — > d<fi/iT. In the range < < tt, is a single-valued function of x. 
Consequently, we can make a change of integration variables in which we replace with x. 
Using Eq. we rewrite x as 



1 



x 



(a 2 + [3 2 - 2a(3 cos 



(B2a) 
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where 



a 



2w 



2z 



(B2b) 
(B2c) 



Note that, owing to the constraint G, a and f3 are real and positive. The Jacobian J for the 
change of the variables is given by 



/ dx 



l 



a 



where 
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(B3b) 
(B3c) 
(B3d) 
(B3e) 
(B3f) 



Using Eq. (jB3a|) . we arrive at the following form for the phase space: 



p poo roc f'\ -z pi 

= de a de b / dy 
J<i> Jo Jo Jo Jo 
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Jo J J 7Ty'C(y-y^)(y + -y) 



(B4) 



where we have interchanged the x and y integrations, and we have not yet specified the 
ranges of integration on the right-hand side of Eq. (jB4|) . 

Now let us work out the ranges of integration. The denominator on the right-hand side 
of Eq. (jB4|) has zeros at y — y+ and y = (There is no zero at C = 0, since y + and ?/_ 



become infinite at that point.) As can be seen from Eq. (|B3a|) . the denominator of Eq. (|B4)1 
is proportional to sin0. Therefore, its zeros correspond to <ft — and <p = n, which are the 
end points of the <ft integration. We conclude that the range of y is restricted to y_ < y < y + . 
In order for y to have a nonzero range, we must have D > 0. This implies that x lies either 
below both zeros of D(x) or above both zeros of D(x). We now argue that the latter range is 
unphysical. Clearly, for x large enough, then Eq. (|B2|) has solutions only for y lying outside 
the physical region < y < 1 — z and Q(y) = 1. Furthermore, y is on the boundary of the 
physical region iff x is equal to a zero of D(x). [To see this, note that the zeros of D(x) 
occur precisely when the range of y vanishes. From Eq. (jB2a|) . we see that the range of y 
vanishes iff a = or f3 = 0. But a = 0or/5 = 0iff?/ = 0or?/=l — 2 or 6 = 0, that is, iff y 
is on the boundary of the physical region.] Then, by continuity, we conclude that y lies in 
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the physical region iff x is restricted to be less than either of the zeros of D. Therefore, we 
require that x satisfy 



x < x max (z,e a ,e 6 ) = Min 



2e a e fe , 



e a + e b 



2z 



(B5) 



Since, in this range of x, y automatically satisfies < y < 1 — z and Q(y) = 1, we 
can drop those explicit constraints on y. For fixed z, the constraint = 1 implies that 
e a + eb > (1 — z)l (2z). This constraint has already been satisfied by virtue of the constraint 
x < x max . However, for purposes of improving numerical- integration efficiency, we can 
explicitly impose it on the range of e&. Taking all of the constraints on the integration 
variables into account, we have for the phase-space integration 
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(B6) 



It turns out that there is no y dependence in the denominator of the integrand and 
that the numerator of the integrand is a polynomial in y of degree two. Therefore, the y 
integration is easily performed: 
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The explicit forms that are needed for the integrand in the present calculation are 
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In carrying out the integrations over z, e a , eb, and x numerically, we made use of the adaptive 
Monte Carlo routine VEGAS (35[. We also checked our results using the built-in numerical- 
integration package in mathematica |26|. Special care must be taken in carrying out the 
x integration numerically. In the limits e a — > or — > 0, the range of the x integration 
vanishes. This can cause a large round-off error in the Monte Carlo integration. We avoided 
this problem by making a further change of variables x = 2e a ebt. The integration interval 
for t is < t < x max /(2e a e fe ). 
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